Chapter 14 Geospatial R - Raster Hydro Analyses
14.1 Introduction
The following activity is available as a template github repository at the following link:
For more:
Load necessary packages and data
#install.packages("whitebox", repos="http://R-Forge.R-project.org")
library(tidyverse)
library(raster)
library(sf)
library(whitebox)
library(tmap)
whitebox::wbt_init()
theme_set(theme_classic())
tmap_mode("plot")# brush <- raster("McDonaldHollowDEM/brushDEMsmall.tif")
#
# wbt_resample("McDonaldHollowDEM/brushDEMsmall.tif",
# "McDonaldHollowDEM/brushDEMsm_5m.tif",
# cell_size = 0.00003925595)
#
# brush <- raster("McDonaldHollowDEM/brushDEMsm_5m.tif")
#
# tm_shape(brush)+
# tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
# tm_scale_bar()tmap_mode("view")## tmap mode set to interactive viewingdem <- raster("McDonaldHollowDEM/brushDEMsm_5m.tif")
dem[dem < 1500] <- NA
tm_shape(dem)+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()## Warning: Currect projection of shape dem unknown. Long lat (epsg 4326)
## coordinates assumed.wbt_hillshade("McDonaldHollowDEM/brushDEMsm_5m.tif",
"McDonaldHollowDEM/brush_hillshade.tif")## [1] "hillshade - Elapsed Time (excluding I/O): 0.12s"hillshade <- raster("McDonaldHollowDEM/brush_hillshade.tif")
tm_shape(hillshade)+
tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
tm_scale_bar()## Warning: Currect projection of shape hillshade unknown. Long lat (epsg 4326)
## coordinates assumed.Prep DEM for hydro analysis
wbt_fill_single_cell_pits("McDonaldHollowDEM/brushDEMsm_5m.tif",
"McDonaldHollowDEM/bmstationdem_filled.tif")## [1] "fill_single_cell_pits - Elapsed Time (excluding I/O): 0.4s"wbt_breach_depressions_least_cost("McDonaldHollowDEM/bmstationdem_filled.tif",
"McDonaldHollowDEM/bmstationdem_filled_breached.tif",
dist = 5,
fill = TRUE)## [1] "breach_depressions_least_cost - Elapsed Time (excluding I/O): 0.52s"filled_breached <- raster("McDonaldHollowDEM/bmstationdem_filled_breached.tif")
## What did this do?
difference <- dem - filled_breached
difference[difference == 0] <- NA
tm_shape(hillshade)+
tm_raster(style = "cont",palette = "-Greys", legend.show = FALSE)+
tm_scale_bar()+
tm_shape(difference)+
tm_raster(style = "cont",legend.show = TRUE)+
tm_scale_bar()## Warning: Currect projection of shape hillshade unknown. Long lat (epsg 4326)
## coordinates assumed.## Warning: Currect projection of shape difference unknown. Long lat (epsg 4326)
## coordinates assumed.## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.wbt_slope("McDonaldHollowDEM/brushDEMsm_5m.tif",
output = "McDonaldHollowDEM/demslope.tif",
units = "degrees")## [1] "slope - Elapsed Time (excluding I/O): 0.5s"slope <- raster("McDonaldHollowDEM/demslope.tif")
tm_shape(slope)+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()## Warning: Currect projection of shape slope unknown. Long lat (epsg 4326)
## coordinates assumed.wbt_d8_flow_accumulation("McDonaldHollowDEM/bmstationdem_filled_breached.tif",
"McDonaldHollowDEM/D8FA.tif",
out_type = "specific contributing area")## [1] "d8_flow_accumulation - Elapsed Time (excluding I/O): 0.24s"d8 <- raster("McDonaldHollowDEM/D8FA.tif")
tm_shape(log(d8))+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()## Warning: Currect projection of shape log(d8) unknown. Long lat (epsg 4326)
## coordinates assumed.## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.wbt_d_inf_flow_accumulation("McDonaldHollowDEM/bmstationdem_filled_breached.tif",
"McDonaldHollowDEM/DinfFA.tif",
out_type = "specific contributing area")## [1] "d_inf_flow_accumulation - Elapsed Time (excluding I/O): 0.37s"dinf <- raster("McDonaldHollowDEM/DinfFA.tif")
tm_shape(log(dinf))+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()## Warning: Currect projection of shape log(dinf) unknown. Long lat (epsg 4326)
## coordinates assumed.## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.wbt_wetness_index(sca = "McDonaldHollowDEM/D8FA.tif",
slope = "McDonaldHollowDEM/demslope.tif",
output = "McDonaldHollowDEM/TWI.tif")## [1] "wetness_index - Elapsed Time (excluding I/O): 0.4s"twi <- raster("McDonaldHollowDEM/TWI.tif")
tm_shape(twi)+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()## Warning: Currect projection of shape twi unknown. Long lat (epsg 4326)
## coordinates assumed.## Variable(s) "NA" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.wbt_downslope_index("McDonaldHollowDEM/brushDEMsm_5m.tif",
"McDonaldHollowDEM/TWId.tif")## [1] "downslope_index - **********************************************************************************"twid <- raster("McDonaldHollowDEM/TWId.tif")
twid[twid > 4000000 | twid <= 0] <- NA
tm_shape(log(twid))+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()## Warning: Currect projection of shape log(twid) unknown. Long lat (epsg 4326)
## coordinates assumed.wbt_d8_flow_accumulation("McDonaldHollowDEM/bmstationdem_filled_breached.tif",
"McDonaldHollowDEM/D8FA.tif")## [1] "d8_flow_accumulation - Elapsed Time (excluding I/O): 0.20s"d_inf <- raster("McDonaldHollowDEM/D8FA.tif")
tm_shape(log(d_inf))+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()## Warning: Currect projection of shape log(d_inf) unknown. Long lat (epsg 4326)
## coordinates assumed.wbt_extract_streams("McDonaldHollowDEM/D8FA.tif",
"McDonaldHollowDEM/raster_streams.tif",
threshold = 6000)## [1] "extract_streams - Elapsed Time (excluding I/O): 0.4s"wbt_d8_pointer("McDonaldHollowDEM/bmstationdem_filled_breached.tif",
"McDonaldHollowDEM/D8pointer.tif")## [1] "d8_pointer - Elapsed Time (excluding I/O): 0.5s"wbt_raster_streams_to_vector("McDonaldHollowDEM/raster_streams.tif",
"McDonaldHollowDEM/D8pointer.tif",
"McDonaldHollowDEM/streams.shp")## [1] "raster_streams_to_vector - Elapsed Time (excluding I/O): 0.4s"streams <- st_read("McDonaldHollowDEM/streams.shp")## Reading layer `streams' from data source `/Volumes/GoogleDrive/My Drive/CLASSES/Hydroinformatics/Hydroinformatics_Bookdown/McDonaldHollowDEM/streams.shp' using driver `ESRI Shapefile'
## Simple feature collection with 37 features and 2 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: -80.49733 ymin: 37.23597 xmax: -80.46388 ymax: 37.25634
## CRS: NAtm_shape(streams)+
tm_lines()+
tm_scale_bar()## Warning: Currect projection of shape streams unknown. Long-lat (WGS84) is
## assumed.points <- read_csv("McDonaldHollowDEM/points.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## lon = col_double(),
## lat = col_double()
## )pointsSP <- SpatialPoints(points, proj4string = CRS("+proj=longlat +datum=WGS84"))
tmap_mode("plot")## tmap mode set to plottingtm_shape(log(twid))+
tm_raster(style = "cont", palette = "PuOr", legend.show = TRUE)+
tm_scale_bar()+
tm_shape(pointsSP)+
tm_dots(col = "red")## Warning: Currect projection of shape log(twid) unknown. Long lat (epsg 4326)
## coordinates assumed.
twidvals <- extract(x = twid, #raster
y = pointsSP, #points
method = "simple")
points$twid <- twidvals
points## # A tibble: 3 x 3
## lon lat twid
## <dbl> <dbl> <dbl>
## 1 -80.5 37.2 89693.
## 2 -80.5 37.2 152567.
## 3 -80.5 37.2 54294.raster_stack <- stack(twi, twid, slope, dem)
raster_values <-extract(x = raster_stack, #raster
y = pointsSP, #points
method = "simple")
points <- cbind(points, raster_values)
points## lon lat twid TWI TWId demslope brushDEMsm_5m
## 1 -80.48355 37.24087 89693.31 -8.086087 89693.31 37.42480 2130.564
## 2 -80.48705 37.24570 152567.28 -9.103719 152567.28 54.68245 2301.624
## 3 -80.48477 37.24697 54293.62 -6.767529 54293.62 30.11527 2430.992